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Abstract 

We compute the radiation reaction force on the orbital motion of compact binaries to the 3.5 
post-Newtonian (3.5PN) approximation, i.e. one PN order beyond the dominant effect. The 
method is based on a direct PN iteration of the near-zone metric and equations of motion of an 
extended isolated system, using appropriate "asymptotically matched" flat-space-time retarded 
potentials. The formalism is subsequently applied to binary systems of point particles, with the 
help of the Hadamard self-field regularisation. Our result is the 3.5PN acceleration term in a 
general harmonic coordinate frame. Restricting the expression to the centre-of-mass frame, we 
find perfect agreement with the result derived in a class of coordinate systems by Iyer and Will 
using the energy and angular momentum balance equations. 
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I. INTRODUCTION 



Since the discovery of the Einstein field equations, gravitational radiation has remained a 
matter of major theoretical interest, and has led to extensive theoretical studies on the nature 
and origin of gravitational wave emission from isolated sources. Approximation methods in 
General Relativity, such as the post-Newtonian (PN) expansion, result in both accurate and 
measurable details of the emission. The prospect of the detection of gravitational waves 
bathing the Earth by ground and space based interferometers LIGO, VIRGO and LISA, 
provides a further impetus to such theoretical investigations. Not only would such detec- 
tors enable an important comparison between astrophysical observations and theoretical 
predictions, but they would also ultimately provide strong tests for General Relativity. 

A favourable potential source for the detectors is the radiation-reaction dominated inspiral 
and eventual coalescence of two compact objects (neutron stars or black holes). For such 
systems, which will undergo hundreds to thousands (depending on the masses) orbital cycles 
in the frequency bandwidth of LIGO and VIRGO, relativistic corrections to the Newtonian 
order in the orbital phasing and wave form play a crucial role in preparing the theoretical 
templates. Indeed, the detection and analysis of these waves in the detectors require at least 
a third post-Newtonian (3PN) correction to both the energy fiux radiated at infinity and in 
the binary's equations of motion. 

The objective of this paper is to compute explicitly the radiation reaction force in the 
equations of motion of a compact binary system at the 3.5PN order ~ (9(c~^) in harmonic 
coordinates (in both a general frame and the centre-of-mass frame). It is well known that 
the leading-order radiation reaction effect occurs at the 2.5PN order ~ 0{c-^) Q i, i, li 
IE IE IE II] (see e.g. Q for a review). We shall, therefore, compute the radiation reaction 
at the IPN relative order, which corresponds, via appropriate balance equations, to the 
IPN corrections in the energy and angular momentum radiated by the system at future null 
infinity (relative to the standard quadrupole formulae). The IPN radiation reaction force is 
also responsible for the dominant effect in the loss of linear momentum, widely referred to 
as the gravitational radiation "recoil". 

Up to the 3.5PN order, the conservative terms in the equations of motion are clearly dis- 
tinct from the non-conservative, radiation reaction, terms. This clean separation manifests 
itself as "even" (Newtonian, IPN, 2PN and 3PN) and "odd" (2.5PN and 3.5PN) orders 
respectively. ^ The 4PN approximation, however, contains both some conservative terms, 
and also, a contribution from the radiation reaction. The former terms are given by some 
"instantaneous" functionals of the source, whilst the latter is associated with the gravita- 
tional wave tails [l^, and is given by a "hereditary" type integral, extending over the past 
history of the source. 

The PN assumption of slow-motion limits the validity of the PN expansion to the so-called 
near-zone of the source (r <C cT, where T is a typical period of variation of the source). An 
important consequence of such a near-zone limitation is that one cannot incorporate directly 
the radiation reaction into the local PN expansion, since the radiation effects depend on the 
boundary conditions imposed on the radiation field at infinity (r ^ a, where a is the size 
of the source), notably the famous no- incoming radiation condition imposed at past null 
infinity, r — > -|-oo with t + r/c = const. At present, two different approaches have been 



^ In the present paper we adopt the terminology that an even (respectively odd) term is one having an even 
(odd) power of 1/c in front. 
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proposed and implemented for treating the problem. 

The first method is based on an asymptotic matching between the PN expansion valid 
in the near zone, and the multipolar expansion for the field outside the source. The match- 
ing occurs in the so-called exterior near-zone of the source, defined by a < r -C cT. The 
asymptotic matching was introduced in this field by Burke and Thorne 0,^]. The radia- 
tion reaction force has been derived by matching up to the 3.5PN order for general matter 
systems j^, and even at the 4PN order, which, as we mentioned above, consists of the 
contribution of tails 0, 0|. The most developed treatment for the exterior multipolar 
expansion is the so-called MPM expansion, which combines the multipolar (M) expansion 
with a post-Minkowskian (PM) scheme The general solution of the matching equation 
between the MPM exterior and PN inner fields has recently been obtained jisl . [iq . In the 
present calculation, we shall parametrise the PN metric by some appropriate "asymptoti- 
cally matched" retarded potentials, which incorporate the 3.5PN radiation reaction effects. 



and are introduced in Refs. [17|, ll8[. These potentials at 3.5PN order result from some 



direct integration of the field equations by means of retarded integrals in the same man- 
ner as in Ref. The end result, which we shall obtain, has already been determined 
by Konigsdorffer, Faye and Schafer [i^ (based on the previous work s EH . li^ ) within the 
framework of the ADM Hamiltonian formalism, and by Pati and Will |23j fsee also 0,11^) 
using their variant iteration of the relaxed Einstein field equations in harmonic coordinates. 

The second method is exclusively applicable to compact binary systems (modelled by 
point particles). It consists of using the known PN expressions for the energy and angular 
momentum radiated at infinity, and of assuming that these fiuxes are balanced by the 
corresponding losses of energy and angular momentum in the binary's local equations of 
motion. This method, based on the e nerg y and angular momentum balance equations, has 
been developed by Iyer and Will l27l |. As shown in Refs. ji^ , the requirement of 
energy and momentum balance determines uniquely the radiation reaction force at 3.5PN 
order in a class of coordinate systems. The residual coordinate freedom is entirely specified 
by two arbitrary gauge parameters at 2.5PN order and by six further ones at 3.5PN. The 
2.5PN parameters assume some specific values in the case of the scalar radiation reaction 
potential of Burke and Thorne at 2.5PN order. They are specified by other values in the 
case of the 2.5PN radiation reaction in harmonic coordinates as calculated by Damour and 
Deruelle ji^li^, |3^. At 3.5PN order, there is also complete agreement, in the sense that 
a unique set of 3.5PN gauge parameters can be determined each time, with the scalar and 
vectorial radiation reaction potentials of Blanchet fl2l], which are valid in some extended 
Burke-Thorne type gauge, with the end result of Pati and Will [23j , who work in harmonic 
coordinates, and with Konigsdorffer et al. jl^, who use ADM coordinates. The method of 
balance equations has also been extended up to 4.5PN order (excluding the tails at 4PN) in 
Ref. iJ. 

The two previous approaches have, therefore, enabled the successful determination of the 
2.5 and 3.5PN radiation reaction terms of the compact binaries' orbital dynamics. Our new 
derivation at 3.5PN order is in complete agreement with the latter works, and, in particular, 
we confirm the earlier result of [23 • The principles of the method followed in Ref. [2^ 
are similar to ours, since both methods are based on the PN iteration of the Einstein field 
equations relaxed by the harmonic coordinate condition. There are, however, important 
differences in the implementation of the asymptotic matching procedure, as well as several 
more minor technical differences. In fact, our own method is justified by the end result of 
the particular matching procedure we use, as given in Ref. Furthermore, though of less 
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relevance in the present problem, the external field in our approach is described by a choice 
of multipole moments which is different from the choice adopted in jisl 2^ 2^ . In addition, 



our treatment of the compact objects model the particles by delta function singularities, 
with the help of the Hadamard self-field regular isat ion, whilst Pati and Will |2j| do not 
implement a regularisation scheme and instead, model each of the bodies by some spherical, 
non-rotating, extended fiuid balls. 

The theoretical framework of the present paper is the 3.5PN equations of motion of a 
general matter system, derived in |18| by a direct FN iteration of the metric in harmonic 
coordinates. As mentioned previously, the metric is expressed as a functional of a particular 
set of non-linear retarded potentials. Such a metric is then specialised to the model of two 
delta function singularities by using the standard prescription for a distributional stress- 
energy tensor in General Relativity. To cure the divergencies associated with each particle's 
infinite self-field, we systematically apply the Hadamard "partie-finie" regularisation [s^, 
S 13) or, more precisely, a specific variant of it defined in [s^ [s^. In addition to the 
divergencies due to the point particles singularities, which are dealt with by Hadamard's 
regularisation, Poisson-like integrals arise at high PN orders and are typically divergent 
at "spatial infinity". However, this problem is technically overcome by the introduction of 
alternative, general solutions to the Poisson equation, in the form of some regularised versions 
of the usual Poisson integral, constructed from a specific finite part procedure called FP. 
This conforms with our definition of the multipole moments in the external field (Ref. 
describes the general formalism which allows one to use such a finite part). Furthermore, 
the recent work [3| shows that the same finite part FP is also to be applied when computing 
the 3.5PN radiation reaction effects. 

When investigating the equations of motion at the 3.5PN level, we will not encounter 
any logarithmic divergencies of integrals similar to the ones found at the 3PN conservative 
level. Such logarithmic divergencies are responsible for an "incompleteness" of Hadamard's 
regularisation in treating the sources' singular nature. Indeed, this incompleteness results in 
the appearance of one physical undetermined parameter, called A, in the 3PN equations of 



motion in harmonic coordinates jl8l. I37l|. or alternatively, the equivalent parameter known 
as the static ambiguity in the 3PN Hamiltonian of the particles in ADM-type coordinates 
[3^ 11^. A recent application of dimensional regularisation, in the framework of which the 
logarithmic divergencies correspond to poles when the spatial dimension d approaches 3, 
yiel ded the numerical value of these parameters and showed that they are indeed equivalent 



j4Ci l41|. ^ A complete calculation of the 3PN equations of motion has also been performed 
using an independent method by Itoh et al. |43l.l43|. Here, for the computation of the terms 
at 3.5PN order, we can use Hadamard's standard regularisation or any of the proposed 
variants of it without encountering such problems. In fact, since there are no logarith- 
mic divergencies at the 3.5PN order, nor any associated poles, the result obtained from 
Hadamard's regularisation is identical to the one arising from dimensional regularisation. 

The plan of this paper is as follows. In Section IHl we present the expressions for the 
metric and equations of motion up to 3.5PN order for general isolated matter systems. 
Section ITTTl considers the specific mathematical model of two compact objects, which will be 
described by point-particle singularities. All the elementary potentials needed for the 3.5PN 
radiation reaction terms are computed in Section HVl Finally, we present our final result for 



^ Other ambiguity parameters, ^, k and ^, present in the radiation field of point particles binaries at 3PN 
order, have also been resolved by means of dimensional regularisation (42i |. 
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the binary's 3.5PN acceleration in Section |3 where we also compare it with the existing 
literature. 



II. FORMALISM FOR GENERAL MATTER SYSTEMS 

This Section presents the 3.5PN equations of motion, which are expressed as a function 
of a particular set of elementary non-linear potentials, and are valid for a general smooth 
hydrodynamical "fluid" system in harmonic coordinates. Thus, we assume (initially) that 
the matter system possesses neither singularities nor black holes, and can be described by 
some Eulerian-type equations involving some high relativistic corrections. 



A. Definition of a set of retarded non-linear potentials 



We begin by stating the result of the direct PN iterative method for the near- zone metric, 
valid for r -C c T, which is parametrised by the retarded potentials (given by some retarded 
integral) introduced in Refs. 17, 3|- Convenience and convention dictate the specific form 
of these retarded potentials. The near-zone metric asymptotically matches to an exterior 
far-zone radiative-type metric in the overlapping exterior near-zone. The exterior metric is 
known from a multipolar post-Minkowskian (MPM) formalism [1^. The matching to the 
PN inner metric yields a solution which is globally defined over a ll sp ace-time (in a formal 
sense of PN expansions) in the harmonic coordinate system 15|, ll6|. The 3.5PN iterated 
metric is given by E"" 
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(2.1a) 
(2.1b) 



(2.1c) 



This metric explicitly involves only some even powers of 1/c. Indeed, the odd terms, no- 
tably the 3.5PN terms we are looking for, are implicitly contained in the definitions of the 
elementary potentials V, Vi, Wij, ■ ■ ■ , which parametrise the metric, as is shown below when 
we perform a PN expansion of the retardation of these potentials. 

The matter stress-energy tensor, T'^'", is conventionally expressed in terms of certain 
mass, current and stress densities, given respectively as. 



J^OO _j_ rj^V 



T 



Oi 



(2.2a) 
(2.2b) 
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(2.2c) 



(where T** = 6ijT'^^). The potentials may be grouped depending on the PN order at which 
they initially appear. For the Newtonian and IPN orders, they are given by, 



V = □^^[-47rGa], 
V, = □^^[-47rG'ai] 



(2.3a) 
(2.3b) 



where D^j^ denotes the usual flat-space-time d'Alembertian retarded (7^) integral. Next, the 
potentials which appear at the 2PN order are defined by. 



X = 



-AtcG Van + W^.d^V + 2VidtdiV + Vd^.V 



Ri — 

-1 



+-{dtVy-2d,V,d,V, 



-AttG {Vai - V,a) - 2dkVdiVk - -dtVdiV 



(2.4a) 

(2.4b) 
(2.4c) 



Finally, the relevant potentials for the highest established order of 3PN are, ^ 
f 



Y — 



-2diVjdjR, - d,V,dtW,, + + 2V,djVd,V + -VAVdy 

+\v'd',V + \v{d,Vf - Imf] , (2.5a) 

-AtxG ( -oRi - aVVi + ^OkWik + ^^^^14 + ^(JkkV^ + WuduVi 
-dtWiAV + d^WkidkVi - dkWudiVk - 2dkVdA - '^VAVd^V 

-^VdtVd,V-2VdkVdkV, + Vd^V, + 2VAdtV, , (2.5b) 

4- = [-AnG V (ai, - 5i,akk) - 2di^^dtV,) + diVkdjVk + dkVid^Vj 

3 

-29(i\4(9fcl^j) - 5ij5feV;„((9fcKn - <9mVfc) - -5ii(<9tF)^ . (2.5c) 



The spatial traces of the potentials will be denoted hyW = Wu and Z = Z^. The harmonic 
gauge condition of the near-zone PN expansion results in four independent PN differential 



When performing the PN iteration for point particles in the context of the extended Hadamard regulari- 
sation js^ilsfil. there are some extra contributions to be added to the 3PN potentials, which are due to the 
violation of the Leibniz rule for the derivative of a product by the distributional derivatives, see Eqs. (3.27) 
in These so-called "Leibniz" terms arise, however, only at 3PN order and do not contribute to the 
present computation at 3.5PN order, where, as previously discussed, the Hadamard regularisation can be 
applied without the problems encountered at the previous 3PN order. 
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identities to be satisfied by the latter potentials, 
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vw + -v"^ 



v, + ^ 



Ri + VV, 



y^ - \%V, + ]^WVi + VRi + y 



R^ + V^V^. 



—Z 
2 



(2.6a) 
(2.6b) 



The harmonicity conditions, Eqs. ()2.6|) . thus provide a verification of certain computed 
potentials which are required at the 3.5PN order. We note, however, that these conditions 
are unable to provide a direct check to the potentials with the most challenging form which 
are necessary for our purposes. 

As suggested by their definitions ()2.3p - ()2.5|) . several recurring structures may be identified 
in our retarded potentials. The advantages of decomposing the potentials are substantial; 
general schemata are developed to solve integrals of a particular form, which not only simplify 
considerably the computational aspect but also reveal interesting analytical solutions. The 
potentials comprise essentially three types of hierarchical terms of increasing complexity 

a) Compact (C) potentials involve spatially compact source terms, proportional to the 
mass, current and stress densities, a, 0", and aij. The support of the source of these 
potentials is limited to the domain of the matter system, which will be given in this 
instance (the black hole "particle-model"), by 5-function singularities (see Section ITTT|l. 
We have, e.g., 

^ \/ = n-i[-47rG'(T], 
= □^M-47rGVaJ, 
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^ijCTfefe)] 



(2.7a) 
(2.7b) 
(2.7c) 



The C potentials are relatively simple, but, for instance, V'^^^ must be calculated with 
the full 3.5PN precision. 

b) Quadratic non-compact (QNC) potentials are generated by spatially non-compact sup- 
ported distribution of the (source-induced) gravitational field. Specifically, such po- 
tentials include terms of the symbolic type ~ n^dVdV , which denote the quadratic 
product of two compact potentials V and/or Vi and their space-time derivatives. Ex- 
amples of QNC potential terms are, 

r{QNC) 



7{QNC) 



^(QNC) 



n^' [-d,vd,v] , 

[-2d^,VdtVj) + d.VkdjVk + dkVidkVj 



3 



-2di^iVkdkVj) - 5ijdkV^{dkV^ - dmVk) - -^5ij{dtVf 

wif^dlv + 2VidtdiV + vdfv 



+-{dtVf -2d,V,d,V 



[2.8a) 



(2.8b) 



(2.8c) 
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c) Cubic non-compact (CNC) potentials include more complicated integral expressions 
of the form of the product of a quadratic QNC potential and a compact C potential. 
The symbolic form reads ~ [^ti (dVdV) dV~\ . The paradigm of such terms is, 

(2.9) 

which must be evaluated at relative 1.5PN order. The other CNC terms, present in 
the potentials T and Yi, will only need to be controlled at the 0.5PN order. 



^(CNC) ^ Q^l 



(QNC);.2 



B. The equations of motion in terms of the potentials 

The 3.5PN equations of motion in the case of a smooth hydrodynamical fluid are obtained 
by replacing the expression of the metric ()2.1|) into the law of covariant conservation of the 
matter stress-energy tensor, VJT^^ = 0, which is equivalent to the equation of geodesies in 
the case of point particle sources. The resulting equation can be expressed as, 

- = n (2.10) 

where P* and P*, introduced here for convenience, can be thought of as some effective linear 
momentum density and force density of the matter system respectively, and are defined by, 

= ^ ^'^""^ (2.11a) 



'Upa j,2 



1 digf^^v'^v" 

2 / iiPt)"' 

'Upa /2 



(2.11b) 



Here, = dx'^/dt denotes the coordinate velocity field {t = x^/c), and we pose v'^ = (c, f*). 
Substituting the metric ()2.1|) into the above expressions ()2.11|) . and performing the FN 
re-expansion, gives 



1 /I 



^2 



+ ^V\' - AVVi + AWiy - 8P, 

c6 Vl6 8 2 ^ 4 

+2#ijt;%2 + 2W^uv%H^ - lOVVy - 20VVjv'v^ 

-Any - 8RjvV + ^Vh' + UVjVjv' + UWijVv^ 
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-mi,V, - 8VRi - + O (^i j , (2.12a) 



F' = diV 



+ ^ (^-V d,V + v"" - 4d,Vj 



i-diV - 2diVj + -VdiV + 2diWj, 
y 8 2 

-AV diVj - 8diRj + ^V^ diV + 8^^- diVj + Ad^X 

+ 4 i—v'^d.V - -9, K v'v^ + —V d,V + diWjk v'^v^v^ 
\IQ 2 ^ 8 

27 

-lOVj diVvh^ - lOVdiVj vh^ - AdiRj vh^ + —V^ diV 



+12VjdiVj v"^ + QWjk diV vh^ + W diWjk vh^^ + Qd,X 
+8diZjk v^v^ - 20VjVdiVv^ - lOV'^ diVj - 8I4. diWjk 
-8Wjk diVk - 8Rj diV - 8V diRj - IQdiYj 

—V"^ diV - AVj Vj diV + 16Rj diVj + 16Vj diRj 

~8V Vj diVj - AXdiV - AV diX + 169^^) + O . (2.12b) 



From this, we deduce the coordinate acceleration as, 



a' 



F' - j^{P' - v') . (2.13) 



The radiation reaction terms in the equations of motion, ()2.12|) - ()2.13|) . appear exphcitly 
when performing the PN expansion of our elementary non-linear potentials. They are ob- 
tained by the careful consideration of all possible contributions at the 2.5PN or 3.5PN orders 
in each of these potentials. More precisely, the 2.5PN and 3.5PN contributions arise both 
from the expansion of the retardation of the inverse d'Alembertian retarded integrals, and 
from the PN corrections already present in the sources of the potentials. For instance, the 
sources, Eqs. (j2.3p - (j2.5|) . involve the mass, current and stress densities, o", 0"^ and aij given 
by (j2.2|) . which depend themselves on the potentials in a manner consistent with the itera- 
tive PN formalism. In addition, some required contributions also occur from the systematic 
order reduction of the accelerations, which is applied when calculating the time derivatives 
associated with the retardations of the potentials, or for instance, when performing the total 
time derivative of the linear momentum density function, in Eq. ()2.13p . By order reduc- 
tion, we refer to the replacement of the acceleration by its explicit expression given by the 
PN equations of motion in terms of the bodies' positions and velocities, and the subsequent 
PN re-expansion, which is performed in a consistent manner at the PN order in question. 

From Eqs. (|2.12p - (j2.13p . it is apparent that the expressions for all the potentials or their 
spatial derivatives are required at a relative 0.5PN order above the existing explicit 3PN 
expansions computed in [l8|. They must be expressed solely as a function of the masses 
and velocities (after the order-reduction of the accelerations). For the moment, Eqs. ()2.12|) - 
()2.13|) were derived in the case of general matter systems, and the precise mathematical 
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description of the source, i.e. the matter stress-energy tensor T^^'^, is required next. Section 
mil introduces the (5-function model of the compact binary and the associated regularisation. 

The PN precision in each of the potentials and their gradients, which we require in order to 
control the 3.5PN acceleration, is given in the following Table. It is convenient to distinguish 
between the computation of a potential and the one of its gradient, because in the case of 
the odd terms, the gradient is often easier to compute. This is due to the simplification that 
the 0.5PN relative term in the expansion of the retardation of a potential (we refer later to 
such an odd-parity term as "retardation-like") is always a mere function of time, which thus 
vanishes when taking the gradient. On the other hand, the gradient of a potential is often 
required at a higher PN order than the potential itself, so it is generally good practice to 
perform a separate computation for the gradient. At 3.5PN order, it is necessary to develop 
the odd terms in. 





to order 


0{c 




diX 


to order 


0{c- 


-3 


V 


id. 


0{c 


-% 


X 


id. 


0{c 


-1 


Vi and djVi 


id. 


0{c- 


-% 


Zij and dkZij 


id. 


0{c- 


-1 


Wij and dkWij 


id. 


olc 


-'), 


Yi and djYi 


id. 


olc 


-1 


Ri and djRi 


id. 


0{c- 


-'), 




id. 


0{c- 


-1 



Apart from the purely compact support potentials, V and Vi, the above potentials consist 
of both compact (C) and non-compact (QNC and/or CNC) support distributed sources. 
The following Sections IIIIHIVI systematically treat how to evaluate each of these different 
types of contributions. We shall find that the evaluation of the required integrals for the 
QNC and CNC terms yield explicit closed-form expressions, valid at any field point over all 
space, for all the odd parts of potentials in the previous Table. This is in contrast with the 
computation of the equations of motion at the previous 3PN order where closed-form 
solutions to certain non-linear Poisson-like integrals could not be given at any field point but 
existed only at the location of each particle (in a regularised sense). As SectionllVlillustrates. 
the alternative "direct" evaluation method of Poisson-like integrals at the location of the 
particles, using the same method as for the 3PN equations of motion [3, |3^ , provides a 
further verification of our analytic closed-form Poisson-like solutions. 



III. APPLICATION TO POINT PARTICLES 

The compact binary system is modelled as two structureless point particles with masses 
mi and m2, which are described by 5-function singularities, and move on a general, not 
necessarily circular, orbit. We neglect the intrinsic rotations (spins) of the particles.^ As 
part of the formalism of fill], we assume that although the equations of motion ()2.10|) - ()2.13|) 
were derived under the assumption of a general smooth stress-energy tensor, i.e. C°°(M^), 
they remain valid in the case of point particles, provided that we supplement the calculation 
by a consistent use of a self-field regularisation. As discussed in the Introduction, the 
Hadamard self-field regularisation is appropriate for the present purpose. 

The use of the point particle model is physically justified in the case of compact objects 
by the fact that, using a Newtonian argument, the tidal effects are formally equivalent to a 



If necessary, the spins may be added to the formahsm along the hnes of Refs. |45|, |46l ]£ 
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correction of the order 5PN ~ 0{c compared to the Newtonian force law (see e.g. (48l|). 



extended compact bodies, as derived in Refs. j23|, l49l ISOL ISlH. a re in complete agreement 
with those derived for the model of point-particles in |l7l. l28l l29l. Is^. The same is true at 
3PN order: there is agreement between the 3PN equations of motion of extended compact 
bodies [i^ 3| and those for point particles [H, |4l[. A general way to justify the use 



of structureless point-particles in order to describe gravitationally condensed objects is to 
invoke the "effacing property" of General Relativity (a consequence of the strong version of 
the equivalence principle). According to this property, the internal structure is effaced when 
considering the motion and the radiation of the compact bodies, so that one can describe 
them only by their masses [30]. Once the use of delta functions is physically justified, the 
advantage, of course, is that they considerably simplify the calculations. 

Following the standard prescription in General Relativity, we write the distributional 
stress-energy tensor of point particles as, 

T'^^(x,t)=/ii(t)<(t)<(t)5(x-yi(t)) + l ^ 2, (3.1) 

where x is the field point, yi and y2 are the positions of the particles, vi(t) = dyi(t)/dt 
is the coordinate velocity in harmonic coordinates, and Vi = (c, vi). The symbol 1^2 
means the same terms but with the particles' labels 1 and 2 exchanged, and 6 denotes the 
usual Dirac three-dimensional delta function. Alternatively, for the mass, current and stress 
densities defined by Eqs. fl2.2|) . we find. 



a = /Ii (5(x - yi) + 1 ^ 2, (3.2a) 
a, = /iiT;i5(x-yi) + l ^ 2, (3.2b) 
aij = fiivivi5{^-yi) + l ^ 2, (3.2c) 

where the quantities /ii and Jli are some explicit functions of coordinate time t, through the 
source trajectories, yi(t) and y2(t), and velocities, vi(t) and V2(t), given by, 

/^iW = / (3-3a) 



v2 



1 + ^ 



(3.3b) 



Here, {gpa)i and {g)i denote the values of the metric and its determinant computed at the 
position of the particle 1 following the prescription of the Hadamard partie-finie introduced 
in Eq. ()3.5|) . As we emphasise later, we do not encounter any problems associated with 
the "non-distributivity" of Hadamard's regularisation at the present 3.5PN order; so, for 
instance, [g gpa)i in Eq. (j3.3|) may be replaced by the product of regularisations, {g)i {gpa)i- 
Thus, it is unnecessary to take into account the subtleties associated with several possible 
choices for the stress-energy tensor in the context of Hadamard's regularisation, which de- 
pend on whether the factors of the delta-function 5(x — yi) are supposed to be evaluated 
at any field point x or at the particle's position yi as assumed in Eqs. ()3.3p . ^ The latter 
problems resulted in the appearance of some ambiguities in the application of Hadamard's 



For instance, a different prescription, valid in the context of the extended Hadamard regularisation, was 
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regularisation at the 3PN order. The ambiguities have since then been resolved by means of 
dimensional regularisation 4^, 41 , i^l , but do not concern us here as we are interested in the 



3.5PN approximation which is unambiguous. We are thus following here the straightforward 
prescription for the stress-energy tensor of point particles in General Relativity. 

We must now be more specific on the method chosen for computing the metric coefficients 
at the point 1 in Eqs. (j3.3|) . Let -F(x) be a typical function we encounter in the problem, 
where for convenience, we indicate only the relevant dependence of the function on the field 
point X (it depends also on coordinate time t through the source positions and velocities). 
The function -F(x) is smooth on M^, except at the singular points yi and y2, around which 
it admits a power-like singular expansions of the type (for any N eN), 

F{^)= E r?/l'^l(n0 + O(rf), (3.4) 

ao<a<N 

(and similarly when 1^2), where ri = |x — yi| — *■ and the coefficients of the various 
powers of ri depend on the unit direction of approach to the singularity, rii = (x — 'yi)/ri. 
The powers of ri are relative integers, a G Z, and bounded from below by some typically 
negative integer oq, depending on the F in question. The coefficients /|"^ for which a < 
are called the singular coefficients of F. The class of functions such as F is called JF. The 
Hadamard partie finie of F G at the singular point 1 is then defined by the angular 
average. 



{F)i = J l^/i Xni), (3.5) 

where dQi = dQ{ni) is the solid angle which is centered on yi and in the direction rii. In 
principle, the Hadamard partie finie is non-distributive (with respect to multiplication) in 
the sense that {FG)i ^ {F)i{G)i in general for F and G belonging to J-'. However, for 
the terms occurring at 3.5PN order, it is unnecessary to account for this feature, as all the 
functions encountered at this order will in fact be such that {FG)i = (F)i(G)i. 

The Hadamard partie finie of an integral, in short Pf J (i^xF, which has divergencies due 
to the singular expansion ()3.4p of the function around at the singular points yi and y2, is 
defined by the always existing limit. 



Pf I d-'xF = lim<! / d^xFix 



^0 



S(s) 

a+3 



a+3<0« + 3 



dQr ff\n,) + In (^^^ J dQ, 



+1^2}, (3.6) 

where in the R.H.S., the integration extends on the domain S{s) =M.^ \ Bi{s) U -62(3), i.e. 
defined by the whole space from which one has excised two coordinate balls Bi{s) and B2{s) 
centred on the two particles and having the (same) radius s. The other terms are defined 



given in Ref. . In this prescription, the factor in Eq. (|3.3a() involving the determinant of the metric is 
calculated at x whilst the other factor is evaluated at yi; in addition, a special version of the Dirac delta- 
function, desi gned in such a way that it permits to keep track of the Lorentz invariance of the formalism, 
is assumed in 1361. 
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with the help of the singular coefficients in the expansion of the function given by ()3.4p . 
The regularisation ()3.6p depends in principle on two constants, Si and S2, appearing in the 
logarithmic terms. These constants play an important role at 3PN order, but will never 
appear in the present work, since there are no logarithmic divergencies at the 3.5PN order. 
In Eq. (j3.6|) . we suppose that the integral converges at infinity, when |x| — > +00. We shall 
explain below how one treats the integral in the case where it diverges at infinity. 

In the extended version of Hadamard's regularisation [i^fiij, one associates to any F E J-' 
a partie finie "pseudo-function", i.e. a linear form defined on the set JF, which permits us 
to give a precise meaning to the notions of Dirac delta functions, and derivatives of singular 
functions in a distributional sense, when they are mutiplied by or act on other singular 
functions in the class J-'. The detailed construction of Ref. jl^ was useful at 3PN order, but 
is not needed in the present work. Nevertheless, it is convenient, because of the availability 
of the computer programmes used in [l^, to adopt all the rules of the extended Hadamard 
regularisation (we know anyway that the different variants of Hadamard's regularisation give 
the same answer at 3.5PN order). In particular, the Dirac delta-function is defined in the 
extended Hadamard regularisation by (for any F G JF), 



Pf / ci3x5(x-yi)F(x) = (F)i 



(3.7) 



where {F)i is the partie finie of the function given by ()3.5p . and where the indication Pf 
reminds us that the equality is true in the sense of the partie finie integral (|3.6p . ^ 

In addition, the derivatives of singular functions, say diF, are to be performed in a 
distributional sense, and for the present work, we use the explicit formula of the extended 
Hadamard regularisation, 

9.F=(9.F)„,,,_^ + D,[F], (3.8) 

where the first term represents the derivative in the ordinary sense (as algebraic computer 
programs would compute it), and where the second term is the purely distributional part of 
the derivative, given explicitly by. 



B,[F]=Ann\ 



^ /hi] 



ni 



+ 5Z •^1 



2-fc] , 



ni 



k>0 



5(x-yi) + l ^ 2. 



(3.9) 



The distributional terms depend only on the singular coefficients of the expansion of F. It 
was shown in Ref. jl^ that this derivative generalises the usual distributional derivative of 
distribution theory in the context of the class of singular functions J-'. It is such that one 
can integrate by parts any integrals; in particular, the integral of the gradient of any F G JF, 
considered in the previous distributional sense, is always zero. Multiple derivatives, as well 
as time derivatives, are treated in a similar way, and the reader is referred to (3^ for details. 
Notice, however, that the distributional derivative ()3.8|) - ()3.9|) . like the usual distributional 
derivative of distribution theory j3^, is seen not to satisfy the Leibniz rule for the derivation 
of a product. This poses a problem at the 3PN order (by the presence of certain ambiguity 
parameters, later resolved by means of dimensional regularisation), but not at the next order 
of 3.5PN. 



^ The "partie finie delta- function" , satisfying H3.7|l . has been developed in Ref. |j35| by the limiting case of 
a particular class of pseudo- functions defined from the notion of the Riesz delta function [s^ . 
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In our investigations, we shall always consider singular functions F G in the form of 
a PN expansion. In order to be systematic, we introduce a special notation for the PN 
coefficients (with odd or even-type parity) of the function F, say, 

F{^,t,c) = J2\Fi^,t). (3.10) 

C (n) 

In the present paper, we neglect the dependence of the PN coefficients on the logarithm of 
c, since such Inc terms do not occur at the 3.5PN order. The coefficients F(n) (and their 
gradients) will be computed at the points 1 or 2 by means of the partie finie ()3.5|) . leading 
to the evaluation of such objects like, 

(n) J 47r („) 

where F stands for any of the PN iterated potentials V ^ Vi, Wij, ■ ■ ■ defined in Section ITTl 



(3.11) 



IV. COMPUTATION OF THE NON-LINEAR POTENTIALS 
A. Compact-support potentials 

Following the nomenclature convention introduced in Section |nj compact (C) support 
potentials refer to inverse d'Alembertian retarded integrals of some source terms, which 
possess as a factor, the source mass, current or stress densities defined by Eqs. fl2.2|) : the 
latter, therefore, for our model of point particles, are of the form Fi 6i and 1^2, where 
Fi(x) G and 5i(x) = 5(x — yi). We shall illustrate the scheme by the derivation of 
the 3.5PN term in the purely compact-support potential V = V^'^^ defined by Eq. (j2.3a|) . 
Performing the expansion of retardations inside the d'Alembertian integral, we obtain, 

^(-'.') - Gg (1)7 d-'xix'- xr-V(x.«) + O (^) . (4.1) 

The source density, a, is substituted by its expression valid for two point-masses, Eq. ()3.2a|) . 
where we recall that Jli{t) is a function of time defined by (|3.3p . This results in, 

(/^i ^1) + 7^76 (/^i ^1) - 1^ (/^i ^1) 



120 c5 dt^ 12Qc'' dt^ 5040 c7 dt'' 

+ 1^2 + 0(1), (4.2) 

where ri = |x — yi| (for convenience we call x the field point in this formula). Notice that 
the 1/c term is a mere function of time, and thus vanishes when taking the spatial gradient 
(hence we employ for this term the notation for a total time derivative d/dt instead of the 
partial derivative d/dt). Furthermore, one can see that this term is actually of order l/c" 
since the Newtonian approximation to /Ii, namely mi, is constant. Performing repeatedly 
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the time derivatives of ri introduces some accelerations which must be consistently order 
reduced by means of the PN equations of motion. In Eq. ()4.2|) . we require the expression 
of Jli{t) up to 3.5PN order, which we easily find by inserting the PN metric ()2.1|) into the 
definition (|3.3|) . and results in, 



mi 





1 




1 


+ 






1 


+ 






1 


+ 





-2 {W)i + \ {V% + \{V)i v\ - 4 mi <^\v\ 



-8 (Z)i - 4 + 2 {W)i - 4 {yi)i (V^,)i - - {y% 
+ Y (V^')i ^i' - 8 (^.)i < + 2(#.,)it;it;i - 3 t;? 

33 1 1 

-4 - 10 v\vl + + Yq < 



(4.3) 



Here, the value of each of the elementary potentials is taken at the singularity 1 following the 
Hadamard partie finie ()3.5|) . For all the terms in Eq. ()4.3|) . the distributivity of Hadamard's 
partie finie at this order is verified, e.g. {WV)i = (W)i {V)i. 

Evidently, the computation of Jli and V proceeds using the PN iteration, where one 
begins with V at Newtonian order, given by [using the notation ()3.10p - ()3.1H) ]. 



V 

(0) 

(0) 



Grrii 

n 

Gm2 
ri2 



1 ^ 2, 



(4.4a) 
(4.4b) 



where ri2 = |yi — y2|- This is then inserted into the IPN term of ()4.3|) in order to obtain 
Jli at IPN order, which hence enables one to deduce V itself at IPN order, and so on. By 
taking into account the fact that there is no odd term in Jli at order 1/c^, we find that the 
first odd term in V at the level is given by. 



V 

(3) 



dt 



(2) 



(0) 



2, 



(4.5) 



(where /i(Q-)i = mi). However, the dominant odd term in the gradient diV is only at order 
1/c^. Notice that there is an odd term ~ 1/c in the case of Wij and for all the potentials 
besides V and Vi. The computation of both {W)i at 0.5PN order and {V)i at 1.5PN order 
are required to get Jli at 2.5PN order, and we have. 



(5) 



mi 



(3) 



2{W)i 

(1) 



(4.6) 
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A useful feature of the potential W = Wu is that it can be expressed at the 0.5PN order 
in a simple way using the following compact-support form [easily deduced from ()2.4c|) ]. 



W = U 



hence we obtain, 



n 



(4.7) 



W 

(1) 



^4 

dt 



d^X \ CTii - - o V 



(0) 



2 (o){o) 



d 



dt 



2Gm^-\vi--{V)i] +1^2 



2 (0)' 



(4.8) 



The 0.5PN term is a purely spatial integral with compact support; it is only a function 
of time, hence (W^(i))i = However, for some more involved potentials, the 0.5PN 

term will also depend on space, because of the contribution from the Poisson integral of a 
corresponding 0.5PN term in the source of the potential. 

Continuing in this manner, and using the explicit computations of the NC potentials as 
explained in the following Sections, we then obtain Jii at the required 3.5PN order. 



(7) 



mi <; ~{v)i - 2(#)i + {VUV)i + \{V)ivl 

(5) (3) (0) (3) 2 (3) 

-A{Vi)iv\ - 8(Z)i - 4(X)i + 2{VUW\ 

(3) (1) (1) (0) (1) 



(1) (1) (1) 



(4.9) 



which allows the compution of V at this order in a straightforward way, 

1 ^ d ,^ ^ Id'^ . 1 93 , 



V 

(7) 



G 



- 1^1 

n (7) 



d , 1 92 ^ 
^1) + o7;t?(^i^iJ 



dt 



(6) 



2dt^ 



(5) 



(4) 



1 55 ^ 



1 

5040 



(4.10) 



(0) 



In fact, we only require the gradient of V at 3.5PN order, and so it is unnecessary to compute 
the second term in the R.H.S. of ()4.10|1 . which vanishes when taking the gradient. All these 
calculations are systematically performed using algebraic computer programmes. 



B. Quadratic non-compact support potentials 

The structure of non-compact (NC) support potentials is that of a d'Alembertian retarded 
integral whose "source" term has a spatially non-compact support distribution. The integral 
is perfectly well defined provided that some sensitive boundary conditions are given for the 
decay of the field at past null infinity (the no incoming radiation condition). However, 
when we expand the retardations of the d'Alembertian integral, some Poisson-like integrals 
will appear at high PN order, which typically become divergent due to the boundary of 
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the integral at (spatial) infinity. This is the well-known problem of divergencies of the PN 
expansion, which is related to the near-zone limitation of the validity of the PN expansion. 

A solution of the latter problem has recently been proposed in Ref. Essentially, the 
work showed that the PN expansion can in fact be iterated ad infinitum by using a 
particular solution of the Poisson equation at each step, which constitutes an appropriate 
generalisation of the usual Poisson integral with a non-compact support source. In this 
specific approach, the source term of the Poisson integral is multiplied by a factor |x|^, 
where B is a complex parameter, and the solution is defined by the "finite part" in the 
Laurent expansion of the Poisson integral when the parameter B tends to zero. In a more 
recent work 16], we have written the end result of [l5| in an alternative form, which shows 
that one proceeds with the expansion of retardations in the PN algorithm by inserting the 
factor |x|^ inside the integrand and taking the finite part in the above sense. This procedure 
will give the correct result for the radiation-reaction odd terms up to the 3.5PN order (we 
are of course still within our specific approach). Starting from the 4PN level, this procedure 
will also have to take into account the appearance of tail contributions in the radiation 
reaction (see Refs. Ullil for the details). 

Hence, we compute the PN expansion of any elementary potential by using such a finite 
part prescription (denoted FPb=o in the following) to cure the problem of divergencies of 
the integrals at the boundary at infinity. Let us consider as an example the computation of 
the non-compact support part of the potential Wij, say, 

i^,f^^)^n^i[-9,v9,v], (4.11) 

which is required, as we have already seen, up to the relative 1.5PN level. By expanding the 
retardations up to this level, we obtain, 

W^r'nx',t) = FP {a- [-r-a.y5,y] +i^A- [-r^d.Vd.V] 

ATTcdt J I I * J 247rc3 dt^ J ' ' ' i * J J 
+0 (1) . (4.12) 

Once again, this particular way of performing the PN expansion, using the regularisation 
FPb=o "at infinity", is justified by our previous works [15,, ^16)]. From this expression, one 
sees that the odd terms at the level 1/c and 1/c^ will come either from the Poisson or 
Poisson-like integrals (which are always even) applied to the odd terms already present in 
the corresponding source, or from the odd terms coming directly from the expansion of the 
retardations (as applied to the even part of the source). We shall henceforth refer to the 
first type of odd-parity terms as Poisson-like, and to the second type as retardation-like; 
"Poisson-like" terms correspond to an even-parity operator applied to an odd source, whilst 



^ We call it the operator of the instantaneous potentials in Refs. [isill^ . denoted by 

+00 / n \ 2k 



fc=0 
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"retardation-like" terms consist of an odd integral (containing explicitly an odd power of 
1/c in front) with an even integrand. 

To begin with, we see that up to 1.5PN order in Eq. ()4.12j) . the odd terms are only 
retardation- like, i.e. given by the two terms with explicit powers of 1/c and 1/c^ (indeed 
the first odd term in the gradient diV occurs only at 2.5PN order). We now compute these 
retardation- like terms; later, at the occasion of more complicated potentials, we shall also 
see how to compute the Poisson-like terms. We insert into Eq. (I4.12|) the precise form for 
the non-compact support source, diVdjV, using the FN expansion ()4.2p up to order 1.5PN, 



c2 



c2 



dj)k {ri) - v\v[d(^i dj)ki (ri) 

t% (^) d,)u in) - v',v% d,),i in) 



2 



+ ^ ( 4 ) • (4.13) 



Since Jli is a mere function of time, it will not affect the subsequent reasonning which deals 
with the spatial integrations. Simply, in the end, we replace Jli by its explicit expression 
at 1.5PN order as deduced from Eq. ()4.3p . Notice again that ()4.13p is purely "even" up to 
2PN order. We have replaced, where appropriate, the masses Jli and JI2 by their Newtonian 
values mi and 

We can distinguish two types of terms in (j4.13|) : as suggested by their names, "self" terms 
refer solely to a single particle and are proportional to e.g. Jli or ml, whereas "interaction" 
terms are functions of both particles and involve for instance the product Jlijl2- Furthermore, 
once the interaction terms are known, the "self" terms can easily be deduced by taking the 
limit y2 —>■ yi (and 1^2). We will, therefore, now focus on how to solve the interaction 
terms (we shall find that the self terms are in fact zero in this calculation). 

In order to compute firstly the term of order 1/c in ()4.12|) to the 1.5PN order, we no- 
tice that each of the partial derivatives in ()4.13p . which acts at the field point x, can be 
transformed into a derivative acting on the source point, either yi or y2, and thus, one can 
merge the factors together and factor the differential operator outside the integral. Hence, 
for instance. 



1 \ . / 1 \ .9 (9 |x 



FP / d^^l^l^di — a- — = —- — r FP / c/^xr^. (4.14) 
B=oJ \rij '\r2j dyldyiB=oJ nr2 

Now, the remaining integral in the R.H.S. of ()4.14|1 can be operationally computed as a 
particular case of the known elementary integral called Yl. In fact, we shall need two 
different types of such integrals, defined by, 

YLiyi,y2) = FP / d^xlxl^^, (4.15a) 

Ztt b=o J rir2 



TL(yi,y2) = -TT Fn / d'^l^^XL-, (4.15b) 

ZTT B=0 J r2 

where L = ii ■ ■ ■ is a multi-spatial index of order i, xl denotes the STF product of spatial 
vectors: xl = STF(x*^ ■ ■ ■ x*'^) also denoted X(^i,) = xl, and the factor — l/27r has been 
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installed for later convenience. These integrals were introduced and computed in [53|, |5J, |5L, 
The general expressions are given in their simplest form by (see 5^ for a detailed derivation), 



Yr 



Tr 



rv2 



E{L-P P) 
yl 2/2 S 



p=0 



12 



3(£ +!)(£ + 2) 



(L-P P) 
2/2 ■ 



(4.16a) 



(4.16b) 



p=0 



Notice that the finite part operator FP b=o plays a crucial role in such calculations; it removes 
any divergency at infinity and makes our computation perfectly clean and well-controlled. 
From the particular case £ = 0, we find Y = ri2 and T = r^2/6j and we thus obtain the 
integrals needed to compute the 1/c term in Eq. ()4.12|) . 

Consider next the retardation-like 1/c^ term in Eq. ()4.12j) : as shown by ()4.12|) and ()4.13p . 
this consists in essence of finding the solution to the integral of the type, 



/ 



d^x Ixl 



Ix'-xP 



rir2 



(i'^x |x 



B 



2x' 



ik 



rir2 



d^y. Ixl 



X 



rir2 



+ X' 



'|2 



\B 



rir2 



(4.17) 



As we see, the second and last terms in the R.H.S. can be computed from Yi with ^ = 1 
and £ = respectively. From ()4.16ap . one obtains Y^ = {y\ + 1/2) f^- The first term, on 
the other hand, involves a priori a new structure, and one defines an elementary function 
for this term, which reads as, ^ 



5L(yi,y2) = -— FP 1x1^+2 

An explicit expression for this elementary function is (see 



XL 

rir2 



(4.18) 



5L(yi,y2) 



ri2 



i£+l){i + 2) 



E(L-P P) 
y\ 2/2 



p=0 



{£+1- p)yl --{p + l){e+l- p)rl, + {p+ l)yl 



(4.19) 



which thus permits us to close our computation of the above retardation-like odd terms. 
(The self terms, obtained from the interaction terms by the limit y2 —>■ yi, are zero.) 

Summarising, we find that our potential W^^^^^ admits the following odd parts. At the 
0.5PN level, it is given by a mere function of time t (through the time dependence of the 
source points yi, y2 and the velocities vi, V2). For later convenience, we introduce a special 
notation for it, and pose. 



W 

(1) 



(QNC) 



(1) 



(4.20a) 



Actually, one easily checks that this function is related to the previous ones by, 
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(4.20b) 



At the next 1.5PN order, the potential is given by a quadratic function of the position, which 
we express into the form. 



W 

(3) 



(QNC) 



(3) (3) (3) 



where the functions of time, as introduced, are respectively given by, 
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(ri2 1/1^ + ^121/2) 



G'^mim2 d^ d'^ ri2 



(3) 



6 cit^ <9i/i*9i//^ ■ 



(4.21) 



(4.22a) 
(4.22b) 
(4.22c) 



As we shall see subsequently, the form of (j4.2ip is of particular use when one derives the 
expression of the crucial cubic-order potential. 



C. Cubic non-compact support potentials 



In order to illustrate the more difficult types of computation which involve both Poisson- 
like and retardation-like spatial integrals, we consider the example of the part of the non 
compact support potential X, referred to as cubic non-compact (CNC) in Eq. ()2.9|) . and 
given by. 



^(CNC) ^ Q^l 



mc)r.2 



(4.23) 



in which we recall that W^f'^'^'^ is defined by fl4.11|) . It is necessary to compute the odd term 



in X*^*-""^'-'^ up to the required order 1.5PN. By expanding the retardations up to this level, 
we obtain, 



FP <! 
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(4.24) 
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The last two exphcit terms in the latter expression represent the retardation-like odd terms 
and are computed using the same techniques as in the case of the quadratic non-compact 
support potential, W-j^^'^^ (see the previous Section). 

In contrast, however, to the case of the QNC potential, we see from Eq. (j4.24j) that since 
the source of the Poisson integral will contain some odd terms at the 1/c and 1/c^ levels, 
there are, in addition to the retardation-like terms, some "Poisson-like" contributions to 
the odd terms in the CNC potential [the first two terms in Eq. ()4.24|) ]. Since there are 
no terms at orders 0.5PN and 1.5PN in the gradient of V, one finds that the 0.5PN and 
1.5PN pieces in the source term, W^f^'^'^df^V, come only from the odd terms in the potential 



W-^j itself, which has already been computed in Eqs. ()4.20|) - ()4.22|) . We show only the 
more difficult case of the 1.5PN order. Hence, the source term for the Poisson integral is in 
the form, 



(QNC) o2 



ai, dl V + aij dl V + Pti df. V+la dl V ■ 



(3) (3f ^^'(0) (l)'" "{2) (3) 



(0) (3) ■ 



(0) 



(4.25) 



We work out this expression using the known even part of the V potential. By thus trans- 
forming the derivatives with respect to the field point into derivatives with respect to the 
source points (using di = —d/dy\ when acting on a function of ri), one arrives at a new 
expression which can be immediately integrated, with the result: 
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(4.26) 



where we recall that the time-dependent coefficients have been given in Eqs. (jOnbjl and 
()4.22|1 . Concerning the other Poisson-like integral, the computation is easier because it is 
required only at 0.5PN order, and we simply obtain. 



FP 

B=0 



B ^mC)Q2y 
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2. 



(4.27) 



Having in hand all the required terms at any field point x, we can now compute the appro- 
priate gradients of these potentials which enter the equations of motion ()2.12|) - ()2.13|) . and 
then obtain their value at the location of the particle 1 using Hadamard's partie finie ()3.5p . 
For instance, the gradient of the term ()4.27p . taken at the location of the particle 1, reads, 

^3 



dk FP A- 

B=0 



B ^(QNC) ^2 



(1) 



ijk 



Gm<2 



• 12 

24 



(4.28) 



where ri2 = |yi — y2| represents the particles' separation. 



D. Alternative derivation of the contribution of NC potentials 

As Sections IIVBI and IIVCI discuss, the previous method consisted of first computing all 
the required NC potentials at any field point x, and then accounting for their contributions 
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to the equations of motion of particle 1 (say), by applying the Hadamard partie finie ()3.5p . 
One can, however, also compute directly the value at point 1 of any potential given in the 
form of a Poisson integral. This alternative technique is indeed used in the computation of 
the equations of motion at 3PN order in Ref. (l^. In the problem of the 3PN equations of 
motion, it was impossible to obtain the value of all the potentials at an arbitrary field point 
X, and thus, the complete result could only be achieved by using the latter technique. In 
the present case of the 3.5PN term, as seen previously, all the potentials at any field point 
could be computed by analytical methods. We are, therefore, able to perform an important 
verification of the result by directly evaluating the potentials "on the line" , i.e. at the source 
point yi (say). 

Such a direct evaluation on the line applies equally well to the computation of retardation- 
like odd terms or the Poisson and Poisson-like integrals. It proceeds from the Hadamard 
partie- finie formalism reviewed in Section 11111 and, in particular, from the expression of the 
partie finie integral (j3.6j) . We check that all the retardation-like odd terms can be cast 
in the form of some sum of partie finie integrals of the type (j3.6p . which are eventually 
multiplied by some spatial vector positions; this is due to the factor, say |x — x'p'^, which 
enters such integrals, and can always be expanded as in Eq. ()4.17|) . The computation using 
this technique of the Poisson-like terms is more complicated but the full proofs are given in 
section|V|of Ref. [3^. Here, we provide only a summary of this computation, which has been 
systematicall performed to obtain full confirmation of our result. Typically, one is dealing 
with the Poisson integral of some F G JF, say, 

P(x') = --^Pf / ^%i^(x), (4.29) 
4tt J |x — x'l 

where the symbol Pf refers to Hadamard's partie finie in the sense of Eq. (13. We are 
interested in the value of this Poisson integral when x' tends to the particle's position yi, 
and this value takes the meaning of Eq. ()3.5|) . In other words, the computation of the 
quantity (-P)i, and also of the corresponding gradient, {diP)i, are required. They have 
already been given in Eqs. (5.5) and (5.17a) of j35|, and read as. 



1 r F 
(P)i = -— Pf / d^^- + 

47r J ri 



In ( ^ ] - 1 

Sl 



[riF),, (4.30a) 

(a,P)i = -l-Fij rf3x^F + ln(^^^ (niriF)i. (4.30b) 

Recall that the partie finie integral ()3.6|) depends on two constants si and S2. In addition, 
r[ = |x' — yi| in ()4.30|) is a "constant" which tends toward zero when evaluating the partie 
finie. Although these constants played an important role in the computation of the equations 
of motion at 3PN order, none of them appear in the present work. For this calculation, we 
also require the formulas concerning the twice-iterated Poisson integral, and the results were 
provided by the equations (5.16) and (5.17b) in [35]. 

We have found by using the above formulas that the direct method is in complete agree- 
ment for all the terms with the analytical methods reviewed in Sections IIV Bl and IIV CI 
However, notice that for the agreement to work, one must crucially take into account, in 
addition to the formulas such as ()4.30p . the contribution of the distributional part of deriva- 
tives. For this calculation, we use formula ()3.9|) . which gives the distributional derivative in 
the extended Hadamard regularisation. 
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V. THE 3.5PN COMPACT BINARY ACCELERATION 



We finally present our result, which gives the complete radiation reaction force in the 
equations of motion of the compact binary at the 3.5PN order, for general orbits and in a 
general harmonic coordinate frame. We write the 3.5PN acceleration of the particle 1, say, 
in the form. 



ai 



+ -A}™ + -Ar + 



2PN 



2.5PN 



3PN 



+ -Af5™ + O 



where the first term is given by the famous Newtonian law, 

Gm2 



ni2. 



(5.1) 



(5.2) 
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The acceleration of the second particle is obtained by exchanging all the labels 1^2. The 
conservative (with even-parity) approximations IPN, 2PN and 3PN have been computed 
elsewhere; they are thoroughly given by Eq. (7.16) in or by Eq. (131) in together 
with the value of the ambiguity parameter A = —1987/3080 computed in |4ll |. 

The result central to our study concerns the radiation reaction (odd-order) acceleration 
coefficients A^^^ and Af^™. We find, for the 2.5PN term,^ 
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(5.3) 



This result is in perfect agreem ent with previous calculations in a general harmonic coordi 
nate frame 
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as well as in the center-of-mass frame [231]. The 
new result is the complete expression of the 3.5PN coefficient in an arbitrary frame which 
we find to be given by. 
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Our notation is ri2 = |yi — y2|, ni2 = (yi — y2)/?'i2, and vi = dyi/dt, ai = dwi/dt for the harmonic- 
coordinates velocity and acceleration (together with \ ^ 2). We pose V12 = vi — V2 for the relative 
velocity. The parenthesis indicate the usual Euclidean scalar product, e.g. (ni2Wi2) = 'n.12 ■ vi2- 
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(5.4) 



Recall that we obtain the results ()5.3|) and ()5.4p directly from ()2.12|) and ()2.13p by summing 
up all the contributions of the regularised values of the potentials. The latter were computed 
following the two (rather independent) methods proposed in Section HVl 

We next give the result for the 2.5PN and 3.5PN relative accelerations in the center-of- 
mass frame. For this calculation, we require the transformation equations for converting 
from the positions and velocities of the particles in the general frame to those in the center- 
of-mass frame. Naturally, the center of mass is defined by the nullity of the binary's mass 
dipole moment at the required PN order. All these relations have been worked out at the 
3PN order in Ref . [H^ . Specifically, the transformation equations take the form, 

yi = [X2 + z/(Xi - X2)V] x + z/(Xi - X2)Qv, (5.5a) 
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For center-of-mass quantities, we denote by x = yi — y2, r = |x|, n = x/r the relative binary's separation 
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y2 = [-Xi + z/(Xi - X2)V] X + z/(Xi - X2)Q V. 



(5.5b) 



The coefficients V and Q are given at 3PN order by Eqs. (3.11)-(3.12) in j56|. For our 
purposes, we require the relative IPN even correction term in V and Q, since we are com- 
puting the 3.5PN radiation reaction force, which is at relative IPN order. In addition, 
we have found that the 2.5PN odd correction term in the latter transformation equations, 
gives a crucial contribution at the 3.5PN order in the radiation reaction force. Inspection of 
Eqs. (3.11)-(3.12) in js^, reveals that the IPN correction term exists only in the coefficient 
V [since Q begins at the level], whilst the odd contribution at order 1/c^ is proportional 
to the velocity, and hence is contained only in Q. Using our previous notation (j3.1(J|) for PN 
coefficients, we have the required terms. 
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We have inserted Eqs. (|5.5|) - (j5.6|) into our general- frame result, in order, therefore, to obtain 
the 3.5PN relative center-of-mass acceleration in the form, 



a = ai — a2 



(5.7) 



All the terms up to 3PN order have already been computed in Eqs. (3.16)-(3.18) of 
The radiation reaction 2.5PN term reads as, 
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(5.8) 



whilst the 3.5PN contribution is given as. 
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(formerly denoted ri2 = r and ni2 = n), and by v = dx/dt = vi — V2 the relative velocity. The mass 
parameters are given as, 

TOl 7712 

m = mi +1712, Xi = — , X2 ^ — , V = XiX2. 

m m 
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4/ 366 _ \ Gm(nvf f 294 376 

12z/ H ^ — u 

V 35 J r V 5 5 

W{nvf (114 + 12z/) - 112{nvy] n} . (5.9) 

The expressions fl5.8|) an d (15.91) allow for an important and mutual consistency check with 
the results of Iyer and Will (26l.l27l|. These authors computed the radiation reaction force at 
3.5PN order for compact binary systems in a class of coordinate systems, by applying the 
energy and angular momentum balance equations to relative IPN order. As a consequence of 
the gauge freedom, the radiation reaction at 2.5PN order depends on two gauge parameters, 
denoted a and /?, whilst at the 3.5PN order, it depends on six further gauge parameters, 
denoted 5a, A = 1, ■ ■ ■ ,6 [l^H^. The eight parameters, a, P, Sa, were proved to correspond 
exactly to the unconstrained degrees of freedom which relate to coordinate transformations. 
For instance, the 3.5PN radiation reaction potentials of Ref. were shown when evaluated 
in the case of com pac t binaries to correspond to a unique, self-consistent choice of all these 
gauge parameters |26l . 127^ . 

In the present case, one also finds perfect agreement between our specific expressions, 
()5.8|) and ()5.9|) , derived from "first principles" , with the end result of Ref. [l^ ■ This 
is provided that the eight parameters, a, P, 6a, assume some constant values reflecting the 
present choice of harmonic coordinates. By comparing our results (|5.8p and (|5.9|) with 
Eqs. (2.12) together with (2.18) in ^2^ (notice the change of notation: £5 — * 6q), we indeed 
obtain a unique and consistent choice for these parameters, given by. 



a = - 1, /5 = 



6, = 'Jl + Qu, S2 = -'i - liy, 

(5.10) 

S3 = fl-f^^, = 10, 

^5 = ^ + ^6 = -t +> 

These values correspond to harmonic coordinates. For these, we find complete agreement 
with the result obtained (also from first principles) by Pati and Will |23|] (see also fH?']). 
On the other hand, some other values correspond to ADM coordinates, as computed by 
Konigsdorffer et al. [20j. Still other values for these parameters, obtained in ji^l, correspond 
to the 3.5PN radiation reaction potentials |l5] valid in extended Burke-Thorne gauge. 

In conclusion, we have computed at 3.5PN order the radiation reaction effect in the local 
equations of motion of a compact binary system in a general harmonic coordinates frame. 
The result was derived using a direct PN iterated expansion of the metric in harmonic 
coordinates at the 3.5PN approximation derived in Refs. 0, 3|- The 3.5PN metric is ex- 



pressed as a function of a particular set of non-linear retarded potentials, defined for any 
general smooth "hydrodynamical" matter distribution. The existence of singular functions 
and divergent integrals, a consequence of our modelisation by two delta-function singular- 
ities, required the Hadamard partie finie regularisation in order to remove each particle's 
infinite self-field. Analytical techniques were used to solve the more complicated integrals 
of the NC supported distribution of the gravitational field, resulting in a final expression 
in closed analytic form. We found that the 3.5PN term in the equations of motion in the 
center-of-mass frame is perfectly consistent with the general expression of the IPN radiation 
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reaction force, derived by energy and angular-momentum balance equations in a class of 
coordinate systems [2^ |23|. This study thus confirms that all the results to date for the 
3.5PN equations of motion, computed either by "direct" iterative PN expansion or by "in- 
direct" energy and angular momentum balance considerations, are fully in agreement with 
each other. Since the equations of motion up to the 3PN order have already been derived 
elsewhere f jisl . 37, 41,0, i^l in harmonic coordinates, j3^ |3^ li^ in ADM coordinates), 
we conclude that the problem of the local equations of motion of compact binaries is solved 
up to 3.5PN order. 
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